
debugDict = {
    'ModSchema': ['QAM','QAM','PSK','PSK'],
    'pulseShaping': ['rrc','square','rrc','square'],
    'txAmp': [40,40,40,40],
    'rxAmp': [20,20,20,20],
    'BER': [0.495,0.59,0.45,0.5545]
}


def plotBaseBandData(data,symbol,rxSymbol=[]):
    colors = np.zeros(len(data))
    if len(rxSymbol)>0:
        figure, axis = ploot.subplots(1, 2)
        for i in range(len(rxSymbol)):
            colors[i] = rxSymbol[i]/4
        axis[1].scatter(data.real,data.imag,c=colors)
        for i in range(len(symbol)):
            colors[i] = symbol[i]/4
        axis[0].scatter(data.real,data.imag,c=colors)
        ploot.show()
    else:
        for i in range(len(symbol)):
            colors[i] = symbol[i]/4
        ploot.scatter(data.real,data.imag,c=colors)
def Convolve(data1,data2):
    ans = np.zeros(len(data1),dtype=type(data1[0]))
    for i in range(0,len(data1)):
        for j in range(0,len(data2)):
            if i-j >=0:
                ans[i] = ans[i] + data1[i-j]*data2[j]
    return ans
def getSNRvsBER(snrMin,snrMax,step):
    snr = np.linspace(snrMin,snrMax,step)
    ber = np.zeros(len(snr))
    for i in range(len(snr)):
        rxAmpDB = 20
        pbRx = np.convolve(pb,ir)
        pbRx = pbRx[0:len(pb)]
        pbRx = comm.awgn(pbRx,snr[i])
        pbRx = pbRx * (10**(rxAmpDB//20))
        bbRx = comm.downconvert(pbRx, int(samplingRate/symbolRate),carrierFreq,samplingRate)
        rxSymbols = comm.demodulate(bbRx,psk)
        ber[i] = comm.ber(comm.sym2bi(rxSymbols,m),comm.sym2bi(symbols,m))
    fig, ax = ploot.subplots(1, 1, figsize=(4, 4), layout='constrained')
    ax.set_xlabel('SNR')
    ax.set_ylabel('BER')
    ax.set_ylim(0,1)
    ax.set_xscale('log')
    ax.plot(snr,ber,'x',ls='-')
    return
def getBeamPattern(dictofDegree):
    ptrn = np.zeros((37,2))
    for i in range(len(ptrn)):
        ptrn[i][0] = -180 + i*10
        try:
            ptrn[i][1] = dictofDegree[-180 + i*10]
        except:
            ptrn[i][1] = 0
    return ptrn
def plotDiffData(rx_range_min,rx_range_max,samples,opt=True):
    rx_range = np.linspace(rx_range_min,rx_range_max,samples)
    berarr = np.zeros(len(rx_range))
    indx = 0
    for i in rx_range:
        env = pm.create_env2d(frequency=27000,tx_depth=60,rx_depth=5,depth=i,rx_range=100,bottom_absorption=0.1)
        arrivals = pm.compute_arrivals(env)
        ir = pm.arrivals_to_impulse_response(arrivals, fs=100000)
        pbRx = np.convolve(pb,abs(ir)[0:int(0.1*100000)])
        pbRx = pbRx[:len(pb)]
        pbRx = pbRx + noise(100000,len(pb),3.6,0.1)
        #pbRx = generateUAN(pbRx,samplingRate,3.6,20)
        #pbRx = comm.awgn(pbRx,20)
        pbRx = pbRx * (10**(rxAmpDB//20))
        bbRx = comm.downconvert(pbRx, int(samplingRate/symbolRate),carrierFreq,samplingRate)
        if opt:
            rxSymbols = comm.demodulate(bbRx,getDemodConst(env))
        else:
            rxSymbols = comm.demodulate(bbRx,psk)
        print("Range: " + str(i))
        berarr[indx] = comm.ber(rxSymbols,symbols)
        print("BER: " + str(berarr[indx]))
        indx = indx + 1
        plotBaseBandData(bbRx,symbols,rxSymbol=rxSymbols)
    print("Average BER: "  + str(np.sum(berarr)/len(berarr)))
def plotDiffData2(rx_range_min,rx_range_max,samples,opt=True):
    rx_range = np.linspace(rx_range_min,rx_range_max,samples)
    berarr = np.zeros(len(rx_range))
    indx = 0
    for i in rx_range:
        env = pm.create_env2d(frequency=27000,tx_depth=60,rx_depth=5,depth=150,rx_range=i,bottom_absorption=0.1)
        arrivals = pm.compute_arrivals(env)
        ir = pm.arrivals_to_impulse_response(arrivals, fs=100000)
        pbRx = np.convolve(pb2,abs(ir)[0:int(0.1*100000)])
        pbRx = pbRx[:len(pb)]
        pbRx = pbRx + noise(100000,len(pb),3.6,0.1)
        #pbRx = generateUAN(pbRx,samplingRate,3.6,20)
        #pbRx = comm.awgn(pbRx,20)
        pbRx = pbRx * (10**(rxAmpDB//20))
        bbRx = comm.downconvert(pbRx, int(samplingRate/symbolRate),carrierFreq,samplingRate)
        if opt:
            rxSymbols = comm.demodulate(bbRx,getDemodConst(env))
        else:
            rxSymbols = comm.demodulate(bbRx,psk)
        print("Range: " + str(i))
        berarr[indx] = comm.ber(rxSymbols,symbols)
        print("BER: " + str(berarr[indx]))
        indx = indx + 1
        plotBaseBandData(bbRx,symbols,rxSymbol=rxSymbols)
    print("Average BER: "  + str(np.sum(berarr)/len(berarr)))
def generateUAN(signal,fs,snr,w):
    #Water noise spectrum
    xf = fftfreq(len(signal),1/fs)
    nff = np.zeros(len(xf))
    nff[1:len(xf)//2] = np.sqrt(50 + 7.5*np.sqrt(w) + 20 * np.log10(xf[1:len(xf)//2]/1000) - 40* np.log10(xf[1:len(xf)//2]/1000 + 0.4))
    nff[len(xf)//2:-1] = np.sqrt(50 + 7.5*np.sqrt(w) + 20 * np.log10(-1*xf[len(xf)//2:-1]/1000) - 40* np.log10(-1*xf[len(xf)//2:-1]/1000 + 0.4))
    nff = np.power(10,nff/20)
    noise_array = np.zeros(len(signal))
    nSignal = comm.awgn(signal,snr)
    noise_array = nSignal - signal
    noise_array_fft = fft(noise_array)
    fnoise_fft = noise_array_fft * nff
    noise = ifft(fnoise_fft).real
    return signal + noise
def noise(sampRate, lenSig, w, s): 
    f = np.linspace(1, int(sampRate/2)-1, int(lenSig))/1e3 # divide by 1000 to convert Hz to kHz
    nt = 17 - 30*np.log(f) # turbulence noise
    ns = 40 + 20*(s - 0.5) + 26*np.log(f) - 60*np.log(f+0.03) #ships noise
    nw = 50 + 7.5*(w**0.5) + 20*np.log(f) - 40*np.log(f+0.04) # wind noise
    nth = -15 + 20*np.log(f) # thermal noise
    tot_noise = nt + ns + nw + nth # total noise in db per Hz
    temp = np.power(10,tot_noise/20)
    phase = 2*np.pi*np.random.randn(lenSig)
    mag = np.sqrt(temp)
    FFT = mag * np.exp(1j*phase)
    noise = np.real(np.fft.ifft(FFT))
    return noise
def getDemodConst(env):
    arrivals = pm.compute_arrivals(env)
    ir = pm.arrivals_to_impulse_response(arrivals, fs=100000)
    ir = ir[0:int(0.1*100000)]
    dataSize = 56*8
    m = 2
    n = 1
    samplingRate = 100000
    carrierFreq = 27000
    bitRate = 80
    symbolRate = bitRate//n
    txAmpDB = 55
    rrcTaps = 6
    #data Modulation
    data =comm.random_data(dataSize, 2)
    symbols = comm.bi2sym(data,m)
    symbols[10:14] = [0,1,0,1]
    symbols[14:18] = [0,1,0,1]
    psk = comm.psk(m)
    bb1 = comm.modulate(symbols,psk)
    pb =  comm.upconvert(bb1, int(samplingRate/symbolRate),carrierFreq,samplingRate)
    #Power Amplification
    pb = pb * (10**(txAmpDB/20))
    
    
    rxAmpDB = 0
    #Received Signal
    pbRx = np.convolve(pb,abs(ir)[0:int(0.1*100000)])
    pbRx = pbRx[:len(pb)]
    pbRx = pbRx + noise(100000,len(pb),3.6,0.1)
    
    pbRx = pbRx * (10**(rxAmpDB/20))

    
    bbRx = comm.downconvert(pbRx, int(samplingRate/symbolRate),carrierFreq,samplingRate)
    ans = np.zeros(2,dtype='complex128')
    ans[0] = (bbRx[10] + bbRx[12] + bbRx[14] + bbRx[16])/4
    ans[1] = (bbRx[11] + bbRx[13] + bbRx[15] + bbRx[17])/4
    print("_______________________________DEMOD CONST__________________________________________________________________")
    plotBaseBandData(bbRx[10:18],[0,1,0,1,0,1,0,1])
    plotBaseBandData(bbRx,symbols)
    print("____________________________________________________________________________________________________________")
    return ans
def demodulate(signal,constellation):
    symbls = np.zeros(len(signal))
    for i in range(len(signal)):
        min_indx = -1
        min_val = 1000
        for j in range(len(constellation)):
            tmp = signal[i] - constellation[j]
            tmp = tmp.real ** 2+ tmp.imag**2
            if min_val > tmp:
                min_indx = j
                min_val = tmp
        symbls[i] = min_indx
    return symbls


def demodulate2():
    bbtmp = np.zeros(len(symbols),dtype='complex128')
    t = np.linspace(0,1/bitRate,samplingRate//bitRate)
    for i in range(0,len(symbols)):
        Q = np.sum(pbRx[i*samplingRate//bitRate:(i+1)*samplingRate//bitRate] * np.sin(2*np.pi*27000*t))
        I = np.sum(pbRx[i*samplingRate//bitRate:(i+1)*samplingRate//bitRate] * np.cos(2*np.pi*27000*t))
        if I > Q:
            bbtmp[i] = psk[0]
        else:
            bbtmp[i] = psk[1]
    return bbtmp


import arlpy.uwapm as pm
import arlpy.comms as comm
import arlpy.plot as plt
import arlpy.signal as signal
import numpy as np
from scipy.fft import fft, fftfreq,ifft
import matplotlib.pyplot as ploot
import math
import matplotlib.axes as axes


pm.models()
rx_range = 1000
tx_depth = 40
rx_depth = 10
print(np.rad2deg(np.arctan((tx_depth-rx_depth)/rx_range)))
env = pm.create_env2d(frequency=27000,tx_depth=60,rx_depth=2,depth=1000,rx_range=100,tx_directionality=getBeamPattern({0:20,10:20,20:20,-20:20,-10:20,30:20,-30:20}),bottom_absorption=0.1)
arrivals = pm.compute_arrivals(env)
rays = pm.compute_eigenrays(env)
pm.plot_rays(rays, env=env, width=900)
ir = pm.arrivals_to_impulse_response(arrivals, fs=100000)
pm.print_env(env)
plt.plot(abs(ir)[0:int(0.1*100000)],fs = 100000)

1.7183580016554572

                name : arlpy
   bottom_absorption : 0.1
      bottom_density : 1600
    bottom_roughness : 0
   bottom_soundspeed : 1600
               depth : 1000
        depth_interp : linear
           frequency : 27000
           max_angle : 80
           min_angle : -80
              nbeams : 0
            rx_depth : 2
            rx_range : 100
          soundspeed : 1500
   soundspeed_interp : spline
             surface : None
      surface_interp : linear
            tx_depth : 60
   tx_directionality : [[-180.    0.]
                        [-170.    0.]
                        [-160.    0.]
                        [-150.    0.]
                        [-140.    0.]
                        [-130.    0.]
                        [-120.    0.]
                        [-110.    0.]
                        [-100.    0.]
                        [ -90.    0.]
                        [ -80.    0.]
                        [ -70.    0.]
                        [ -60.    0.]
                        [ -50.    0.]
                        [ -40.    0.]
                        [ -30.   20.]
                        [ -20.   20.]
                        [ -10.   20.]
                        [   0.   20.]
                        [  10.   20.]
                        [  20.   20.]
                        [  30.   20.]
                        [  40.    0.]
                        [  50.    0.]
                        [  60.    0.]
                        [  70.    0.]
                        [  80.    0.]
                        [  90.    0.]
                        [ 100.    0.]
                        [ 110.    0.]
                        [ 120.    0.]
                        [ 130.    0.]
                        [ 140.    0.]
                        [ 150.    0.]
                        [ 160.    0.]
                        [ 170.    0.]
                        [ 180.    0.]]
                type : 2D

psk

array([ 1.+0.j, -1.+0.j])


print(len(ir))

# Number of samples in normalized_tone
N = len(abs(ir)[0:int(0.1*100000)])

yf = fft(ir[0:int(0.1*100000)])
xf = fftfreq(N, 1 / 100000)

plt.plot(xf, np.abs(yf))

139


#data Generation and Config

dataSize = 2000
m = 2
n = 1
samplingRate = 100000
carrierFreq = 27000
bitRate = 80
symbolRate = bitRate//n
txAmpDB = 20

rrcTaps = 6
#data Modulation

symbols =comm.random_data(dataSize, 2)
#symbols = comm.bi2sym(data,m)
psk = comm.psk(m)
bb1 = comm.modulate(symbols,psk)
pb =  comm.upconvert(bb1, int(samplingRate/symbolRate),carrierFreq,samplingRate)
pb = pb * (10**(34/20))

#Power Amplification

pb = pb * (10**((txAmpDB)/20))
pb = pb + 1.127
#Plotting

time = signal.time(pb,100000)
plt.plot(time,pb)


pb2 = np.zeros(len(pb))
bb2 =  comm.upconvert(bb1, int(samplingRate/symbolRate),0,samplingRate)
pb2 = bb2.real * np.cos(2*np.pi*27000*time) + bb2.imag * np.sin(2*np.pi*27000*time)
pb2 = pb2 * (pb[0]/pb2[0])


time = signal.time(pb,100000)
plt.plot(time,pb2)


rxAmpDB = 0
#Received Signal
pbRx = np.convolve(pb2,abs(ir)[0:int(0.1*100000)])
pbRx = pbRx[:len(pb)]
pbRx = pbRx + noise(100000,len(pb),3.6,0.1)
#pbRx = generateUAN(pbRx,samplingRate,3.6,20)
#pbRx = comm.awgn(pbRx,20)
#Power Amplification

pbRx = pbRx * (10**(rxAmpDB//20))

#plotting
time = signal.time(pbRx,100000)
plt.plot(time,pbRx)


#Demodulating Signal
bbRx = comm.downconvert(pbRx, int(samplingRate/symbolRate),carrierFreq,samplingRate)
rxSymbols = demodulate(bbRx,getDemodConst(env))
comm.ber(comm.sym2bi(rxSymbols,m),comm.sym2bi(symbols,m))

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________

0.0


plotBaseBandData(bbRx,symbols,rxSymbols)


#High Depth Cases
plotDiffData(150,1500,30)

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 150.0
BER: 0.0

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 196.55172413793105
BER: 0.0

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 243.10344827586206
BER: 0.0

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 289.6551724137931
BER: 0.0

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 336.2068965517241
BER: 0.0

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 382.7586206896552
BER: 0.0

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 429.31034482758616
BER: 0.0

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 475.8620689655172
BER: 0.0

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 522.4137931034483
BER: 0.0

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 568.9655172413793
BER: 0.0


#Without Training, Shallow Water
plotDiffData2(50,150,8,opt=False)

Range: 50.0
BER: 0.0

Range: 64.28571428571429
BER: 0.0

Range: 78.57142857142857
BER: 0.0

Range: 92.85714285714286
BER: 0.001

Range: 107.14285714285714
BER: 0.0

Range: 121.42857142857143
BER: 0.0

Range: 135.71428571428572
BER: 0.1705

Range: 150.0
BER: 0.0

Average BER: 0.0214375


#With Training, Shallow Water
plotDiffData2(50,150,8)

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 50.0
BER: 0.0

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 64.28571428571429
BER: 0.0

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 78.57142857142857
BER: 0.0

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 92.85714285714286
BER: 0.0

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 107.14285714285714
BER: 0.0

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 121.42857142857143
BER: 0.0

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 135.71428571428572
BER: 0.0

_______________________________DEMOD CONST__________________________________________________________________
____________________________________________________________________________________________________________
Range: 150.0
BER: 0.0

Average BER: 0.0


!jupyter nbconvert --to html equalizer.ipynb

[NbConvertApp] Converting notebook equalizer.ipynb to html
[NbConvertApp] Writing 5121214 bytes to equalizer.html